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The mathematical model of a steadily propagating Saffman-Taylor finger in a Hele-Shaw channel 
has applications to two-dimensional interacting streamer discharges which are aligned in a periodic 
array. In the streamer context, the relevant regularisation on the interface is not provided by surface 
tension, but instead has been postulated to involve a mechanism equivalent to kinetic undercooling, 
which acts to penalise high velocities and prevent blow-up of the unregularised solution. Previous 
asymptotic results for the Hele-Shaw finger problem with kinetic undercooling suggest that for a 
given value of the kinetic undercooling parameter, there is a discrete set of possible finger shapes, 
each analytic at the nose and occupying a different fraction of the channel width. In the limit in 
which the kinetic undercooling parameter vanishes, the fraction for each family approaches 1/2, 
suggesting that this ‘selection’ of 1/2 by kinetic undercooling is qualitatively similar to the well- 
known analogue with surface tension. We treat the numerical problem of computing these Saffman- 
Taylor fingers with kinetic undercooling, which turns out to be more subtle than the analogue 
with surface tension, since kinetic undercooling permits finger shapes which are corner-free but not 
analytic. We provide numerical evidence for the selection mechanism by setting up a problem with 
both kinetic undercooling and surface tension, and numerically taking the limit that the surface 
tension vanishes. 

PACS numbers: 47.15.gp 47.20.Ma 52.80.Mg 
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I. INTRODUCTION 

Inspired by the seminal work of Saffman and Taylor 
[I] , an enormous amount of research has been undertaken 
on aspects of the problem of a steadily moving finger of 
inviscid fluid in a Hele-Shaw cell of channel geometry 
(for an overview of Hele-Shaw flows, with a thorough 
discussion on flows in the channel geometry, see Hi]). 
In a typical experiment mm, air is injected from the 
left end of a horizontal channel, which is otherwise filled 
with viscous fluid. The air-fluid interface is unstable, 
as the less viscous fluid is displacing the more viscous 
fluid (the Saffman-Taylor instability). As the interface 
evolves from left to right, a fingering pattern develops, 
which ultimately results in a single finger of air propa¬ 
gating steadily along the Hele-Shaw cell and occupying a 
fraction A G (1/2,1) of the channel width. As the finger 
speed increases (via higher injection rates), the ratio A is 
observed to decrease towards roughly A = 1/2 [1]. 

The most common approach to understanding the 
structure of the corresponding mathematical model is 
to study the role of a dimensionless surface tension pa¬ 
rameter cr, which decreases as the finger speed increases 
[6]. There are exact solutions for the special case tr = 0 


[SISIIZI, but these all take the ratio A as an input param¬ 
eter, and so do not describe the observed experimental 
behaviour. The ‘selection’ of A = 1/2 as the physically 
appropriate solution in the limit ct —>■ 0 is a difficult prob¬ 
lem in exponential asymptotics [EHUj. The relevant anal¬ 
ysis predicts that, for a given tr, there is a discrete set 
of solutions with I/2<A<1. Ascr—>-0, the value of 
A for each solution branch approaches the special value 
A = 1/2. Numerical solutions support these conclusions 

[51 ITS]. 

In the present study, we are concerned with the effect 
that kinetic undercooling has on the Hele-Shaw problem 
in a channel geometry. The appropriate dimensionless 
model for a steadily propagating finger is [16] 



= 0 

in 

^00 ■) 

(la) 

d(j) 

dn 

= 0 

on 

5Hoo, 

(lb) 

</> = 

X 

= i_A 

on 


(Ic) 

d(j) 

dy 

= 0 

on 

y = ± 1 , 

(Id) 

X 




(fr- 

1-A 

as 

X -T —00, 




A < I 2 /I < 1, 

(le) 

(fr- 

—X 

as 

X -T -boo. 




- 

1 < y < 1. 

(If) 


scott.mccue@qut.edu.au 





2 


Here (j) is the velocity potential in the frame of reference 
of the finger, djdn denotes a directional derivative nor¬ 
mal to the interface dilao, Vn is the normal velocity of 
the interface, and c is the kinetic undercooling parame¬ 
ter. The unregularised version (zero kinetic undercool¬ 
ing, c = 0) of Eqs. Q has Eq. (Ic I replaced by 


1-A 


on 


( 2 ) 


which also applies for the zero surface tension case men¬ 
tioned above. 

Kinetic undercooling-type conditions arise in a variety 
of applications. In the Hele-Shaw context, the kinetic 
undercooling term arises from the curvature in the trans¬ 
verse direction (perpendicular to the parallel walls of the 
Hele-Shaw cell), and its dependence on the interface ve¬ 
locity. This effect was included by Romero who 
modelled the contact angle as a linear function of the ve¬ 
locity, leading to a boundary condition such as Eq. ( [T^ . 
An alternative interpretation is to consider the existence 
of a wetting layer of the receding fluid that remains on 
the plates of the Hele-Shaw cell. Park and Homsy [18] 
derived a power-law relationship between the thickness 
of this layer and the capillary number. This relationship 
leads to a power-law dependence on velocity, with the 
term cvn in Eq. (Ic I replaced by cz;)), where 7 = 2/3 is the 


exponent derived in HHj. Such a term may be referred to 
as representing nonlinear kinetic undercooling. The theo¬ 
retical short-time existence of solutions to Hele-Shaw flow 
with this regularisation was established by Pleshchinskii 
and Reissig m- Recently, the stability of an expanding 
circular bubble with both surface tension and nonlinear 
kinetic undercooling has been considered, in both linear 
[ 20 IET], and weakly nonlinear |22l [23] regimes. In this 
paper, however, we consider linear kinetic undercooling 
( 7 = 1 ) only. 

In the context of melting or freezing, Stefan-type for¬ 
mulations may include a Gibbs-Thomson law with ki¬ 
netic undercooling [24II27j , with much attention given to 
instabilities and pattern formation at the interface of a 
growing dendrite |5SJ UHl HH] ; in that case, in the limit of 
vanishingly small specific heat, the governing equations 
reduce to those for Hele-Shaw flow. Thus the unstable 
Hele-Shaw model describes the manner in which a super¬ 
cooled liquid freezes, with Eqs. Q above relevant for a 
single dendrite propagating with constant velocity in a 
channel. Kinetic undercooling conditions also apply on 
interfaces in very similar moving boundary problems de¬ 
scribing mass transfer situations, such as the diffusion of 
solvent through glassy polymers [SHI IM] . 

Of particular interest here, the model Q has appli¬ 
cations to streamers, which is a topic that has received 
much attention in the physics literature in recent times 
(see the review [31]). Streamers are finger-shaped elec¬ 
trical discharges which occur during the early stages of 
electric breakdown in sparks or lightning, for example. 
They are caused by subjecting a weakly ionized gas to 
a strong electric field, leading to an ionization reaction 


via collisions of highly energetic electrons with neutral 
molecules. The streamers themselves are characterized 
by a thin charge layer and associated ionization front 
that forms the finger shape. 

A minimal model for streamer discharges consists of 
a coupled system of reaction diffusion equations for the 
electron and ion density. A further equation relates the 
Laplacian of the electrostatic potential 4> to these densi¬ 
ties. For negative streamers, these equations can be ap¬ 
proximated by a moving boundary problem by assuming 
the ionization layer is a sharp interface that separates 
the strongly ionized streamers from the weakly ionized 
gas ahead of front. The result is Laplace’s equation for 
the electrostatic potential outside the interface. For the 
case in which there is a periodic array of two-dimensional 
streamers with equal spacing, all propagating in the x- 
direction with a constant electric field E = —xi in the 
far field as a: —^ 00 , one can impose Neumann conditions 
to isolate a single streamer [3^34] . Under this periodic 
geometry, if the electric field or periodic spacing is suf¬ 
ficiently small (strong interaction between neighbouring 
streamers), the streamers evolve from their initial con¬ 
ditions to a travelling wave profiles, so that they propa¬ 
gate uniformly. The approximate model is then given by 
Eqs. 0. 


In the context of streamers, the boundary condition ® 
has been used instead of Eq. (Ic) (see [33], for exampl^. 
The former is appropriate if the streamer is assumed to 
be ideally conducting (</> = 0 in the streamer) and the 
electric potential is assumed to be continuous across the 
interface. Indeed, the condition ([^ was used by Luque 
et al [33] in their study of periodic streamers (see also 
Ref. [31] )■ However, as is known from the Hele-Shaw 
literature, the unregularised time-dependent model is ill- 
posed, with a dense subset of all initial conditions leading 
to finite time blow-up that is characterised by infinitely 
sharp cusps on the interface [7]. Such behaviour is not 
physical (in either the Hele-Shaw or streamer context). 
The regularising term (Ic) is postulated by Ebert and 


coworkers [36II38| for streamers, and used, for example, 
to model perturbed translating circles [3S1I1S]- A further 
relevant discussion is contained in Ref. [33] ■ Here the 
kinetic undercooling parameter c is proposed to account 
for the thickness of the ionization front. In the present 
paper, we shall employ the language of Hele-Shaw flows, 
but keep in mind the application of streamers, discussing 
the relevance of the analysis and results in Sec. |IV| 


The Saffman-Taylor problem with kinetic undercool¬ 
ing, described by Eqs. 0 , has received modest atten¬ 
tion compared to the surface tension analogue mentioned 
above. The selection problem was treated by Chap¬ 
man and King [16] , who used exponential asymptotics to 
show that discrete families of analytic fingers exist, with 
the finger width for each family tending to 1/2 in the 
limit that the kinetic undercooling parameter c vanishes. 
These authors showed that A ~ 1/2-1- as c —?► 0 for 

each branch, but did not compute the constant a. More 
recently, a numerical study by Dallaston and McCue [JT] 
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showed that, for a given kinetic undercooling parameter 
c, a continuous family of corner-free finger solutions exist 
with widths A G [Amin, !)■ Further, it was found that the 
minimum width Amin —>■ 0 as c —>■ 0. While this continu¬ 
ous spectrum of solutions appears to be at odds with the 
analysis of Chapman and King IE], the two studies need 
not contradict each other since the numerical scheme in 
Ref. [H] is not designed to distinguish between solutions 
with analytic fingers and those with fingers that are also 
corner-free but may not be analytic (that is, for all the 
numerical solutions found in Ref. mi, the first derivative 
exists at the nose, but the higher order derivatives may 
not exist there). 

In this paper we aim to reconcile these results by 
constructing numerical solutions to Eqs. Q that have 
analytic fingers. The rigorous results of Tanveer and 
Xie [42l|43| suggest that solutions to the Hele-Shaw prob¬ 
lem with sufficiently small values of the surface tension 
coefficient must have interfaces that are analytic. With 
this in mind, our strategy is add surface tension to the 
model Q , so that Eq. (Ic I is replaced by 

X 

(l) = aK + CVn- - -r On i 9 floo, ( 3 ) 

i — A 

where a is the surface tension coefficient and 7 is the 
curvature of the interface (see Ref. [U] for an in-depth 
study of Hele-Shaw flows with surface tension and ki¬ 
netic undercooling). Our hypothesis is that the work of 
Tanveer and Xie carries over to Eqs. (la)-(lb), (Idl- 
(§ so that solutions to the problem with kinetic 
undercooling and surface tension must be analytic at the 
nose. Thus with kinetic undercooling fixed at some value 
c > 0 , by taking the limit tr —)■ 0 , we select the analytic 
solutions studied in Chapman and King [Ej. Using this 
strategy, we are able to produce a plot of finger widths A 
versus kinetic undercooling c for the first two branches, 
thus filling in the gap left by Chapman and King |16j and 
Dallaston and McCue m- Our results have implications 
for the problem of periodic streamers studied by Tuque 

et al [55] . 

Our numerical scheme is based on a boundary integral 
formulation, as outlined in Sec. |llj Sec. Ill 


summarises 


our main results, while Sec. |IV| includes a discussion. 


II. BOUNDARY INTEGRAL FORMULATION 

For the formulation of the problem, we follow the 
work of McLean and Saffman | 6 ] and Chapman and King 
0111 ]. Since ((> is a harmonic function, we define an an¬ 
alytic complex potential w{z) = y) -|- ?/), where 

is a, stream function and z = x + iy. The conformal 
transformation z w maps the fluid region onto an infi¬ 
nite strip of unit width in the potential plane. A second 
conformal map, w,—>'X = ^-|-i ?7 = e”’^™ maps this strip 
onto the upper half y-plane. The interface is mapped 
onto the unit interval on the real line, 0 < < 1 , with 

the upper wall mapped onto —00 < ^ < 0 and the center 
line y = 0 mapped onto 1 < < 00. 


The complex velocity can be written 


dw 

dz 


= qe 


-w 


(4) 


where q is the velocity tangential to streamlines, and 0 is 
the angle the tangent to the streamlines make with the x- 
axis. The logarithm of this velocity, log q — i9, is analytic 
in the upper half X'Pl&ne, and its real and imaginary 
parts can be related by a property of Hilbert transforms 
called the Kramers-Kronig relations, such that 


log§ 


1 m - TT 

Wo e - C 




0 < e < 1, (5) 


since 6 = tt everywhere on the real line except the unit 
interval. Note that the integral is of Cauchy principal 
value type. 

Relating the quantities q{^) and 9{^) to the curvature 
of the interface (Ref. 0) allows us to rewrite the dynamic 
condition (§ as the differential equation 

(1 - X)q =(1 - A) 7 r 2 CT( 7 C^ 

-I- cnq^ cos 9— — cos0, 

dC 

We now have Eqs. ^ and ^ relating 
associated boundary conditions 

5(0) = ^, 0(1) = ^, §(1) = 0, (7) 

which correspond to uniform flow at the tail (^ = 0 ) and 
a stagnation point at the nose (^ = 1 ). 

Given values of the physical parameters a and c, we 
seek to solve Eqs and then compute the finger 

width A via 


0 < C < 1. (6) 

q and 0 , with the 


log(l-A) 


1 9{C) - TT 

^ Jo e 


d^'. 


( 8 ) 


which comes from setting ^ = 0 into Eqn (§. 

We now introduce another variable substitution that 
simplifies the equations and removes the explicit depen¬ 
dence on A. We let 9{^) = 9{^) — tt, q{^) = (1 — A)$(^) 
and introduce new parameters 


(T7r^ CTT 

T^’ 2(1-A) 


(9) 


where 7 is a scaled surface tension parameter 00 and 
e is a scaled kinetic undercooling parameter [16] . Then 
with some manipulation, the governing equations become 




.d 6 » 


logq = -- 

TT 


■Jo 


de 

0(0 


d9 


q = iq^^ { f 4" + cos 0 


de 


-C) 


d^', 


( 10 ) 

( 11 ) 
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both of which hold on 0 < ^ < 1 , together with boundary 
conditions 

0 ( 0 ) = 0 , g( 0 ) = l, 0 ( 1 ) = -|, <?(!)= 0 . ( 12 ) 


Given a solution for 0(^), we can calculate the width of 
the finger using Eq. ([^, which is now 

log(l - A) = - r ® dC', (13) 

TT Jo ? 

and calculate the shape of the interface from 


x{0 + ^yi0 = 


1 - ^ exp (i 0 (CO) , 


(14) 


Using McLean and Saffman’s exact solutions for the 
unregularised problem [ 6 ], 


q = 



9 = cos ^ q, 


(15) 


where a = (2A—1)/(1 —A)^ is arbitrary, and the formulae 
for the physical coordinates implicit in Eq. (141, we can 


recover the analytic formula for the shape of the finger 
given by Saffman and Taylor, namely that x(^) = ((1 — 
X)/tt) log^, y{^) = ( 2 A/ 7 r) cos“^ v^. Combining the two 
results gives 


X = 


2(1-A) 


'“«('“(i!))' <“> 


which is often referred to in the literature as the ZST 
solution, being equivalent to the expression derived first 
by Zhuravlev |45j and then by Saffman and Taylor [T] 
(see [5] for an alternative derivation). 



(a) 


A 



(b) 


FIG. 1: (Color online) Numerical results for kinetic 
undercooling e = 0.1. (a) The shape of the fingers for 
the surface tension values 7 = 0.03, 0.5,1 (from 
innermost to outermost curve), (b) The dependence of 
the finger width A on surface tension 7 . Numerically 
computed data points are indicated by the solid (blue) 
circles. The open (red) circle is the estimated finger 
width for this family at 7 = 0 . 


III. NUMERICAL RESULTS 


We solve our system of integro-differential equations 


( 10 )-( 12 | by applying the numerical scheme outlined in 


the Appendix. The approach involves dividing the do¬ 
main 0 < ^ < 1 into N + 1 unevenly spaced grid points 
and solving a system of A — 1 equations for the unknown 
function 9 at each of the N — 1 interior points using a 
Newton solver. The other quantities of interest can be 
computed subsequently. 


A consequence of discretising the integral in Eq. (11) is 


that the A ^—1 equations depend on the unknown function 
9 at all of the grid points, which leads to a fully dense 
Jacobian J in the Newton scheme. In order to proceed 
with a large number of grid points, we have employed a 
Jacobian-free Newton-Krylov method which does not re¬ 
quire the formation of the full Jacobian; instead, a sparse 
approximation is all that is required for preconditioning 
of the Krylov subspace linear solver, as described in the 
Appendix. 

Typically, for a fixed surface tension parameter 7 > 0 
and kinetic undercooling parameter e > 0 , the scheme 


converged to a solution that corresponds to a particu¬ 
lar finger shape with a single finger width A. The initial 
guess used for Newton’s method was either the exact so¬ 
lution (15) for 7 = 0, e = 0, or an already converged 
solution with similar parameter values. For moderate to 
large values of 7 , A^ = 3000 grid points were used, while 
for small values of 7 we used a larger number of grid 
points, up to a maximum of A^ = 5000. 


Some representative finger shapes are presented in 
Fig. |l(a)| Here we have fixed the kinetic undercooling 
parameter to be e = 0.1 and provided results for three 
different surface tension values, 7 = 0.03,0.5 and 1. We 
observe that the fingers are qualitatively the same in each 
case, and that the finger width A is greater than 1/2 and 
decreases as the surface tension 7 decreases. 

Each of these three solutio ns co rrespond to a single 
data point on the curve in Fig. 1 1(b)] which shows the de¬ 
pendence of the finger width A on the surface tension 7 for 
e = 0.1. This figure clearly demonstrates the trend that 
as surface tension decreases, the finger width decreases. 
For values of surface tension below roughly 7 « 0.015, we 
were unable to compute sufficiently well converged solu- 
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tions (using up to iV = 5000 grid points). The reason for 
this breakdown in the numerical scheme is related to the 
singular nature of the limit 7 —>■ 0 , which is illustrated 
by the highest derivative in Eq. (10 1 being multiplied by 
7 . We discuss this issue further below. 

Recall that our hypothesis is that all these fingers are 
analytic curves, since non-zero surface tension does not 
allow non-analytic solutions. On the other hand, for 7 = 
0 (zero surface tension), Dallaston and McCue [H] show 
that there is a continuous family of corner free solutions 
for A > Aniin(e)j where for e = 0.1 the minimum value is 
roughly Amin ~ 0.44. To select a single solution in this 
family (with e = 0.1 and 7 = 0 ) that has an analytic 
finger, we propose to consider the branch of solutions for 
e = 0.1 and 7 > 0 and take the limit 7 —)■ 0 +. 

Since it is difficult to calculate solutions for extremely 
small values of 7 , we use an extrapolation approach to 
obtain an estimate for the finger width at 7 = 0. One 
option to achieve this is to fit a polynomial to the last 
few data points and extract the value of this polynomial 
at 7 = 0. However, we have the result in the case of 
zero kinetic undercooling that A ^ | / 37 ^'^^ as 7 —)■ 0+, 

thus it seems reasonable to suggest that the same scaling 
holds in the case of finite kinetic undercooling. As such, 
we use the relation 


A ~ a-I-,07^/^, (17) 

and fit a small number of the final few points to this 
equation. The value obtained for a is the predicted finger 
width for 7 = 0 , the intercept on the vertical axis in 
Fig. |l(b)[ 

In addition to the branch of solutions shown in 
Fig. |l(b)[ we have found evidence of additional solution 
branches. This is precisely the same behaviour as known 
to occur for the case without kinetic undercooling (e = 0 ) 
[a [min]. Romero m and Vanden-Broeck |15j demon¬ 
strated the existence of multiple solution branches for a 
given value of 7 numerically, and Chapman [S] and others 
proved the existence of an infinite number of branches us¬ 
ing exponential asymptotics. Kessler and Levine [47l l48] 
suggested that only the lower branch is stable while the 
other, higher branches are unstable mm- 

Thus for this particular example e = 0.1, we postu¬ 
late there are a countably infinite number of solutions 
branches, each more difficult to compute than the previ¬ 
ous. We show three such curves in Fig. [2(b)[ Each follows 
the trend of decreasing A as 7 decreases. It is difficult to 
compute A values for small values of 7, but again, we are 
able to extrapolate to estimate the analytic solution for 
7 = 0 on a second branch. For the third branch, the 
lowest 7 value at which the numerical scheme converged 
was too large to give an accurate extrapolation estimate. 

Also shown in Figs. ia) and [^c) are three solutions 
branches for e = 0 and e = 0.2, respectively. Of course, 
the e = 0 case is the original surface tension problem 
[giHiiii]. The extrapolation technique was used on two 
branches for each of these values to obtain an estimate 
for an analytic solution 7 = 0 . In principal we could 


A 



(a) e = 0 

A 



(b) 6 = 0.1 

A 



FIG. 2: (Color online) Dependence of finger width A on 
surface tension 7 for fixed values of kinetic 
undercooling. In each case, three solutions branches are 
shown. The open (red) circles represent an extrapolated 
value for 7 = 0 . 
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construct a similar figure for any fixed value of kinetic 
undercooling, e. 

To provide further insight into the singular nature of 
the limit 7 —>■ 0, we have presented in Fig. plots of the 2- 
norm condition number, cond (J), of the Jacobian versus 
the surface tension 7 for the cases presented in Fig. 

On this log-log scale, the data appears to be linear as 
7 —>■ 0, which implies that cond(J) ~ const 7 “^, where 
p is a positive constant depending only on the kinetic 
undercooling parameter e and the particular branch of 
solution. These observations support the claim that the 
problem is singular in the limit 7—^0 for e > 0, which is 
consistent with our hypothesis that the findings of Tan- 
veer and Xie [HI [H] do not extend to the case 7 = 0, 
e > 0. The singular nature of the problem also helps to 
explain the numerical findings m of a continuous solu¬ 
tion space for 7 = 0 , e > 0 . 

By extrapolating our numerical results for the lower 
two solution branches for many different values of e, we 
have constructed the data provided in Fig. |4(a)| These 
are our estimates of the finger widths associated with the 
analytic solutions to Eqs. also analysed in Chapman 
& King m using asymptotic techniqu es. It is notewor¬ 
thy that our solution branches in Fig. 4(a) also appear 
to approach A = 1/2 in the limit e —)■ O"'', which agrees 
with Chapman & King. 


Also included in Fig. |4(a)| as a dashed curve is the lower 
bound of all solutions, including non-analytic fingers, as 
found by Dallaston and McCue [41]. We see that as e 
increases, this lower bound appears to asymptote to the 
lower solution branch for analytic fingers. 


In Fig. |4(b)| we include more details of the primary 
branch, showing the dependence of A against 0 < c < 1. 
Recall that as e —>■ 00 , c —>■ 1“ and A —1“. Since our 
method is most useful for investigating the primary few 
branches in the (A, e) solution space, we shall not attempt 
to match our curves to the solution curves in m, which 
are only valid near e = 0 in the limit that eN —>■ 0, where 
N is the branch number, that is when A — 1/2 ~ 0{1). 
This implies that their results are only valid for the higher 
order branches. Unfortunately, it is therefore infeasible 
to use our proposed method to investigate these solution 
branches. 


While our main focus is selection as e —>■ O’*', there are 
interesting results in the limit that the kinetic undercool¬ 
ing parameter e —> 00 , or equivalently, as c —> 1 “, which 
we can use to test our approach. Chapman and King 
determined in the appendix of m that 1 — A <C 1 — c 
and that the asymptotic behaviour of the first branch is 
given by 


c ^ 1 — (1 — A) log (1/(1 — A)), as A —1 . (18) 

See the inset in Fig. |4(b)] for a comparison of the numer¬ 
ical results with this asymptotic relation. The shape of 
the finger in this limit is given by Chapman & King [16j 
as being circular at the nose. See Fig. for a comparison 
between this asymptotic solution and solution profiles for 


cond(J) 



(a) e = 0 


cond(J) 



cond(J) 



FIG. 3: (Color online) The condition number of the 
Jacobian for solutions computed with N = 1000 nodes, 
plotted against 7 for the first branch (circles), second 
branch (diamonds) and third branch (squares), for fixed 
values of e. The (blue) dashed line is a rough linear fit 
through the data on the first branch for small e. 




















7 


A y 
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(b) 


FIG. 4: (Color online) This figure constitutes our main 
result, (a) The selection of two distinct branches shown 
as the solid (red) circles. Note the distinction between 
these and the continuous solution space found by 
Dallaston and McCue [41], bounded below by Amin(e), 
the dashed (black) curve. The primary branch seems to 
asymptote to the curve Amiii(e) as e —)■ oo. 

(b) The entire primary branch from part (a). Recall 
that as e —>■ oo, c —>■ 1“ and A —>■ 1“. The portion of the 
branch shown in part (a) is boxed in the lower left 
corner for reference. The inset shows a comparison with 
Eq. (18) provided by Chapman & King Hi, shown as 
the smooth (blue) curve. 


small surface tension and varying values of kinetic under¬ 
cooling. 

We end this section by mentioning the results pub¬ 
lished in Ref. [li. The model Q was treated in 
Ref. |49j and numerical results were presented with dis¬ 
crete branches of solutions. However, as discussed in de¬ 
tail by Dallaston and McCue jH] , these discrete branches 
were due to very small numerical errors, which were cor¬ 
rected in Ref [IT] . 


IV. DISCUSSION AND CONCLUSIONS 

We have computed multiple discrete analytic solution 
branches for the Saffman-Taylor finger with kinetic un¬ 
dercooling, corresponding to those predicted asymptoti¬ 


FIG. 5: (Color online) The shape of the finger for 
solutions on the primary branch with 7 = 0.03 and 
e = 0,1,10, 3500 (from innermost to outermost curve). 
For comparison, the dashed (red) curve shows a 
semi-circle of unit radius. 


cally by Chapman and King |16] . The greatest numerical 
challenge is to distinguish analytic solutions from non- 
analytic ones, given the inability of a finite difference 
scheme to capture high derivatives. Here we achieved 
this goal by extending the numerical scheme of 
to include both surface tension and kinetic undercooling, 
and then extrapolating to find the limit as surface ten¬ 
sion goes to zero for fixed kinetic undercooling values. 
Our numerical results agree with asymptotic results in 
selecting a finger width of 1/2 as e —>■ O+j as well as in 
producing a semi-circular interface as the finger width 
tends to the channel width. 

The inclusion of surface tension ensures that the nu¬ 
merical solutions we compute represent analytic fingers; 
the theory of Tanveer and Xie [42] for the pure sur¬ 
face tension problem ensures that any solution that is 
(which finite differences can certainly distinguish) is 
also analytic, and it reasonable to assume this carries 
over when both surface tension and kinetic undercooling 
are present. The inability of the numerical method to 
distinguish discrete solution branches in the absence of 
surface tension (as observed in Ref. m) suggests that 
the results of Tanveer and Xie do not apply when 7 = 0; 
that is, there do exist but nonanalytic travelling finger 
solutions for the pure kinetic undercooling problem. 

We have not considered the numerical computation of 
the time-dependent version of Eqs. Q. Analytic travel¬ 
ling finger solutions are only relevant if analyticity is pre¬ 
served in evolving from an initial condition. While this 
is unlikely to occur for sufficiently large kinetic under¬ 
cooling (Dallaston and McCue jH] have numerical and 
asymptotic evidence of corner formation for c > 1 ) it 
may be possible if kinetic undercooling is small enough 
(c < 1). Extrapolating a time dependent solution with 
zero surface tension and nonzero kinetic undercooling 
from one with nonzero surface tension and nonzero ki¬ 
netic undercooling may introduce further complications 
given the structural instability of the time-dependent 






















problem in the zero surface tension limit gim]. Any 
numerical scheme would have to be very precise, but also 
avoid the node-crowding effect typical of numerical con¬ 
formal mapping methods. 

We close with remarks about the relevance of our re¬ 
sults for the study of streamer discharges. For this appli¬ 
cation, it has been proposed that Hele-Shaw type mod¬ 
els can be used to approximate the dynamic evolution 
of streamers, with a kinetic undercooling term used as 
a form of regularisation, where the kinetic undercool¬ 
ing parameter is a measure of the actual thickness of 
the ionization front [32l 136144(1] . Recall that Luque et 
al. [33] considered a periodic array of strongly interacting 
streamers and showed that, after some transient period, 
they propagate uniformly. By isolating a single trans¬ 
lating streamer, they treated the Hele-Shaw problem (fl) 
as an approximate model, except that they used Eq. ) 
instead of ([^. That is, they considered the unregularised 
version of the classical Saffman-Taylor finger problem [1] . 
Here we have treated Eqs. 0 with nonzero kinetic under¬ 
cooling, and presented results that support the hypothe¬ 
sis that the width of each streamer finger for vanishingly 
small kinetic undercooling (vanishingly small thickness 
of the ionization front) is one half the period of the array 
of periodic streamers m- This conclusion explains why 
the exact solution to the unregularised problem with the 
free parameter A set to 1/2 agrees with time-dependent 
solutions to the full streamer problem, at least near the 
tip of the streamer [33] • 

As our study suggests, the use of a kinetic undercooling 
type regularisation for evolving streamers is not without 
complications. While the Hele-Shaw model without reg¬ 
ularisation is ill-posed, and therefore not appropriate for 
streamer discharges (or any application, for that matter), 
the time-dependent version of Eqs. Q is still difficult 
to handle numerically. For example, the time-dependent 
version of Q is highly unstable; linear stability shows all 
modes of perturbation (of a flat interface) are unstable 
[HUSO]- Further, it would presumably require a partic¬ 
ularly sophisticated numerical scheme to distinguish be¬ 
tween time-dependent solutions with analytic fingers and 
those that are non-analytic but corner-free. As such, it 
seems that a better dynamic model for streamers may 
involve kinetic undercooling plus another regularisation 
effect that comes from the full streamer model. This ad¬ 
ditional effect may then act like surface tension does in 
the Hele-Shaw context described here, allowing for selec¬ 
tion of physically appropriate solutions to the streamer 
problem of interest. 
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Appendix: Numerical scheme 


We seek to solve the integro-differential set of equa¬ 
tions ( |Io| - 0 and associated boundary conditions ( |12[ ) 
numerically, in a manner similar to McLean and Saffman 
[6| (see also Refs [H] [HI |49|). In order to deal effec¬ 
tively with the integral in Eqn. (11), we note that both 9 
and q are non-differential functions of ^ at the endpoints, 
with square root type singularities at ^ = 0 (the tail) and 
^ = 1 (the nose). The variable transformation 


r = 1 - c' 


(A.l) 


is used to ensure that both variables have at least two 
derivatives at the end points, and 0<r<l/2is the real 
root of the transcendental equation 

-|- 2eT = cot ttt, (A. 2) 


which is obtained from considerations regarding the dif¬ 
ferentiability of 9 at both endpoints [^. 

The integral in Eq. (11) is a Cauchy Principal value 
integral; we can add and subtract the singular part to 
give 


log? = 




0(0 


e 

In 


0 ( 0 - 0(0 




i-e 


(A.3) 


Using Eq. (A.l), the first integral in Eq. (A.3) becomes 

A C'9{C) 


2 


1-C'- 


dC'. 


Since now 9 = 0 at = 1, the integrand has a 
removable singularity there, and ca n be replaced by 
— (1/2) d0/dC|^^;^. Again using Eq. (A.l), the second 
integral in Eq. (A.3) becomes 


C' 


oia-0{o 


(1 - C'2)1-1A (1 - ^'2)l/r _ (1 _ ^2)1/, 


dC', 


which has a removable singularity at = (. L’Hopital’s 
rule is again used to replace the integrand at (/' = ^ with 

-(r/2)d0/dC. 

Now turning to the numerical scheme itself, we dis- 
cretise the unit interval ( S [0,1] using TV -|- 1 nodes, 
0 < Cn ^ 1 where n = 0,1,2,..., A, and look to solve for 
the vector of unknowns u = [^i, 02, ■ • • j Given 

an initial guess Uq for the values of 0„, or an updated 
vector Ufc, we can calculate the values [?i, ? 2 , • ■ •, ? 7 V-i]^ 
using Eq. (11) (rewritten in terms of (/), then substitute 
both 9 and q into Eq. (10) using third order mixed finite 
difference formulas to approximate the derivatives. Thus 
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FIG. 6 : The structure of a typical Jacobian, for 
N = 100. The plot shows log^g |J|, for 7 = 0.02, e = 0.2, 
on the first solution branch. 


we have a system of iV — 1 nonlinear algebraic equations 
for the iV — 1 unknowns in u, which we solve using a 
Jacobian-free Newton-Krylov method m, implemented 
by the SUNDIALS package KINSOL [^. Once con¬ 
verged, the solution can be used to recalculate q, X and 
the physical coordinates x and y. 


Due to the global nature of the integral equation ( |11[ ) , 
the Jacobian matrix J of the nonlinear system is fully 
dense. The Jacobian-free Newton-Krylov method avoids 
the need to form this dense matrix, leading to consid¬ 
erable efficiency gains. It does so by using a precondi¬ 
tioned Krylov subspace solver at the linear level, which 
requires only an approximation of the true Jacobian for 
preconditioning purposes. To efficiently construct this 
approximation, we observe that the largest entries in 
the Jacobian matrix are contained within a narrow band 
around the main diagonal - a consequence of the finite 
difference approximation of the derivatives in Eq. (lOl; 
other relatively large values are located in the rightmost 
columns. An example of this striking pattern is provided 
in Fig. where we see the magnitude of the entries 
in J decay with distance from the main diagonal. To 
construct the preconditioner, we retain only the entries 
within the narrow band and a relatively small number 
of the rightmost columns, yielding a sparse approxima¬ 
tion that is efficient to form and factorise. By varying 
the bandwidth, the trade-off between the cost of factori¬ 
sation and the effectiveness of the preconditioner can be 
controlled. This approach is analogous to that applied 
recently by Pethiyagoda et al. |2S] , who also solved a 
coupled system of two integro-differential equations, de¬ 
rived using a boundary integral method. Similar tactics 
for constructing sparse preconditioners from dense Jaco- 
bians have been implemented for other non-local systems 
(see Ref. |56]). 
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